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We have studied the numerical solutions of Richardson equations of the BCS model in the limit of 
large number of energy levels at half-filling, and compare them with the analytic results derived by 
Gaudin and Richardson, which in turn leads to the standard BCS solution. We focus on the location 
and density of the roots, the eigenvalues of the conserved quantities, and the scaling properties of 
the total energy for the equally spaced and the two-level models. 
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I) Introduction 

One of the most studied models in Condensed Matter Physics is the BCS model of superconductivity, proposed by 
Bardccn, Cooper, and Schrieffer in 1957 pL The BCS model is usually formulated in the grand canonical ensemble, 
which is appropiated to systems with a macroscopic number of fermions. However for small systems, such as nuclei or 
ultrasmall metallic grains, one has to consider the canonical ensemble where the BCS wave function is not adequate. 

Fortunately enough, there is a simplified version of the BCS model which is exactly solvable in the canonical 
ensemble. This occurs when the strengh of the pairing interaction between all the energy levels is the same. The 
exact solution of the reduced BCS model was obtained by Richardson in 1963 [§, who also studied its consequences 
and properties in a series of papers in the 60's and 70's || 0, |j. The exactly solvable BCS model, in turn, is closely 
related to the rational spin model of Gaudin, defined in terms of a set of commuting Hamiltonians [Tc|| . 

The integrability of the reduced BCS Hamiltonian was proved in 1997 by Cambiaggio, Rivas, and Saraceno (CRS) 
Llf, who constructed a set of conserved quantities in involution, commuting with the BCS Hamiltonian and closely 

More recently Gaudin's trigonometric model have been generalized, 
There are also bosonic pairing Hamiltonians satisfying the same 



related to the rational Gaudin's Hamiltonians 
to include a g-term a la BCS, in references [12 
type of equations as in the fermionic case (?|, 14 
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15 1 . Generalizations of the SU(2) BCS model to other Lie groups |16|] 
and supergroups |17[ have been worked out. Both, the BCS and Gaudin's models are intimately linked to conformal 
field theory |L8], |19|, p0|, Chern-Simons theory Jig , and Integrable Models in Statistical Mechanics [|l], [2^, ^3| (see 
ref. for a short review on these topics). The exact solution has also been used to study the effect of level 
statistics j2f| and finite temperature |^6|, |2?J on the ultrasmall metallic grains. 

An important property of the exact solution is that the energy of the states and the occupation numbers agree, 
to leading order in the number of particles, with the BCS theory jjj. This result was obtained by Gaudin jlo| and 
Richardson ||, using an electrostatic analogy of the Richardson's equations, which have to be solved for finding the 
eigenstates of the BCS Hamiltonian. More recent applications of this electrostatic model to nuclear pairing can be 
found in ref. g§. 

In the limit of large number of levels N, Richardson's equations become integral equations, which can be solved 
with techniques of complex analysis JlO| ]. This is sufficient to determine the location and density of roots, as well as 
the total energy of the ground state to leading order in N. To study higher order corrections of the ground state 
energy and the excitations, Richardson developed a multipole expansion of the fields appearing in the electrostatic 
analogue model || . We shall follow Gaudin's approach which is more geometrical, but use also some of the results 
obtained by Richardson. 

The aim of this paper is to make a revision and precise comparison between the numerical solution and the analytic 
expressions for the location and density of roots of the Richardson's equations, as proposed originally by Gaudin in 
his unpublished paper in 1968, which got finally published in his collected papers in 1995 |L0]]. The comparison will 
be made in the limit of large N for the two-level model and the equally spaced model. 

These models have been extensively used in Nuclear Physics since they were introduced by Hogaasen-Feldman p9[ | 
(two-level) and Richardson himse lf |p| (e qually spaced). The latter model has also been used to describe the physics 
of ultrasmall metallic grains ]3(], ftl[ 32|, |33| Q (for a review see ref. p5[). In the large N limit the equally spaced 
model is superconducting for all values of atractive BCS couplings g, while the two- level model displays a quantum 
phase transition between a superconducting state and a normal state as a function of the BCS coupling constant. 

The organization of the paper is as follows: in sections II and III we give a primer of the Richardson's exact solution 
of the BCS model, together with its integrability. In section IV we introduce the electrostatic analogue model of BCS. 
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In section V we review Gaudin's continuum limit of the Richardson's equations and their general solution, using 
complex analysis. In section VI we derive the gap and chemical potential BCS equations under the assumption that 
the roots of Richardson equations collapse into a single arc, and study in detail the two-level and equally spaced 
models. In section VII we present our numerical results and compare them with the analytic ones obtained in section 
VI. We also study the scaling properties of the roots and the ground state energy. In section VIII we present our 
conclusions and prospects. Finally, we describe in an appendix the conformal mappings associated to the equally 
spaced model. 



II) Richardson's exact solution of the reduced BCS model 

The reduced BCS model is defined by the Hamiltonian [Q, [30[ [}5| 

Hbcs = 2 E £ 3<? c ]* c 3° - G ^Z c )+ c )- c 3'- c i'+ > W 

j,v=± j-J' 

where Cj t ± (resp. cj ± ) is an electron annihilation (resp. creation) operator in the time-reversed states \j, ±) with 
energies £j/2, and G is the BCS dimensionful coupling constant. The sums in (|l|) run over a set of N doubly 
degenerate energy levels Ej/2 (j = 1, . . . , N). We adopt Gaudin's notation, according to which e.j denotes the energy 
of a pair occupying the level j We shall assume that the energy levels are all distinct, i.e. e% ^ £j for i ^= j. The 
Hamiltonian (lj) is a simplified version of the reduced BCS Hamiltonian where all couplings have been set equal to a 
single one, namely G. Hereafter we shall refer to (|l|) simply as the BCS Hamiltonian. 

Richardson had long ago solved this model exactly for an arbitrary set of levels, Ej [§, ||, 0- To simplify matters, 
we shall assume that there are not singly occupied electronic levels. As can be seen from (|l]), these levels decouple 
from the rest of the system; they are said to be blocked, contributing only with their energy £j/2 to the total energy 
E. The above simplification implies that every energy level j is either empty (i.e. |vac)), or occupied by a pair of 
electrons (i.e. cj + cj _ |vac)). Denote the total number of electrons pairs by M. Then of course M < N. The most 
studied case in the literature corresponds to the half-filled situation, where the number of electrons, N e — 2M, equals 
the number of levels N |3jJ. In the absence of interaction (i.e. G = 0), all the pairs occupy the lowest energy levels, 
forming a Fermi sea. The pairing interaction promotes the pairs to higher energies and even for small values of M or 



G the levels are pair correlated |3U 32, B3| 



In order to describe Richardson's solution one defines the hard-core boson operators 

h <•, <:,.■■ &} = cj !+ ct_, jy,- = b]b j , (2) 

which satisfy the commutation relations, 

[bj,b],] = Sj, f Q.-2Nj). (3) 
The Hamiltonian ([!]) can then be written as 

"<>, ■< X M> >'?y • (4) 

Richardson showed that the eigenstates of this Hamiltonian with M pairs have the (unnormalizcd) product form 

M N 

\M) = [] B*|vac>, B V =Y, —JT b J < ( 5 ) 

V—l J—l J 

where the parameters E v (y = 1, . . . , M) are, in general, complex solutions of the M coupled algebraic equations 



N 



M 

E 



Ei. — E,, 



(6) 
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which play the role of Bethe ansatz equations for this problem |fj| [L9|, ^C], El], |2^, ^3|, |24j . The energy of these states 
is given by the sum of the auxiliary parameters E v , i.e. 

M 

E{M) = Y,E V . (7) 

The ground state of Hbcs is given by the solution of eqs. (||) which gives the lowest value of E(M). The (normalized) 
states (||) can also be written as H 

\M) = -7= E V'(ii J --- J jM)6t i ---6t M |vac), (8) 

where the sum excludes double occupancy of pair states and the wave function ip takes the form 

M 

The sum in (p|) runs over all the permutations, V, of 1, • • • , M. The constant C in (^) guarantees the normalization 
of the state If (i.e. (M\M) = 1). 

Ill) Integrability of the reduced BCS Hamiltonian 

A well known fact about the BCS Hamiltonian is that it is equivalent to that of a XY spin model with long range 
couplings, and a "position dependent" magnetic field proportional to £j. To see this, let us represent the hard-core 
boson operators (||) in terms of the Pauli matrices as follows, 

bj=cr+, b] = aj, AT,- = 1(1-^), (10) 
in which case the Hamiltonian (Ji|) becomes 

Hbcs = H XY + ~ £ e 3 + G(N/2 - M), (11) 

3 

Hxv = -Y,^-^(T + T-+T-T+), 

3 

where the matrices 

N 1 
T a = ]T^, (0 = 0,+,-), with t d = -a], tf = af, tj = aj , (12) 

3=1 

satisfy the SU{2) algebra, 

[T a ,T b ]=f?T c , with f+° = f°r=-l, f+~=2, (13) 

whose Casimir is given by 

T • T = T°T° + -(T + T~ + T~T + ). (14) 
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This spin representation of the BCS model is the appropiate one to study its integrability, as it was shown by 
Cambiaggio, Rivas, and Saraceno (CRS), who constructed a set of operators, Rj (j = 1, . . . , N) [ pi , 

N t t 

R i = -t<l-2gY J h -^ L i (i = l,...,n), (15) 

satisfying the properties 

[H B cs,Ri] = [Ri,Rj] = 0, (i,j = l,...,N). (16) 
Let us denote by Xj the eigenvalue of Rj acting on the state i.e., 

R j \M) R = X j \M) R . (17) 

CRS, unaware of Richardson's solution, did not give an expression of Xj in their work. This was found in reference 
|0| using CFT techniques, 



[ v i=i(^) 



which verifies J2i — {M — N/2) upon using eqs. so that it vanishes at half filling. 

However CRS showed that Hxy given in eq. (|ll]) can be expressed in terms of the operators Ri as 



3 



Then using the eqs. (|TT|) , (|l8[), and (|[) one can show that E(M) is indeed given by eq. @j. (For the two-level system 
the last term in @ must be replaced by -2G(N/4 + l)JV/4 ). 



IV) The Gaudin-Richardson's electrostatic model of BCS 



The equations (|6j) admit a 2 dimensional electrostatic analogy, which will be very useful in the study of the limit 
N — > oo |To). Let us consider a set of N charges —1/2 fixed at the positions Sj on the real axis, and a uniform 
field parallel to this axis with stregth — 1/2G. The problem is to find the equilibrium positions of M charges +1 at 
positions subject to their mutual repulsion, the attraction with the —1/2 charges, and the action of the uniform 
field. The 2D electrostatic potential is given by W + W* , where W reads, 



W({EJ, {£,}) = \ ]T log( £i \og{E v - Efj) 



2 



(20) 



It is easy to see that the Richardson's eqs. (^) arise as the stationary conditions dW/dE^ — 0, while the conserved 
quantities (|l8|) are proportional to the forces exerted on the fixed charges, i.e. 



dW 

Xi = 2G (21) 

The total energy (Q) is the center of gravity of the charges E^, which must be located symmetrically with respect to 
the real axis, i.e. if E^ is a solution to the eqs. (0), then E* must also be a solution. This condition is fullfilled in two 
cases: i) either E^ is real, or ii) E^ and E* form a complex conjugate pair. The formation of these complex pairs is 
usually an indication of the superconducting properties of the ground state. 
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V) Continuum limit of Richardson's equations 

In reference Gaudin proposed a continuum limit of the equations (|^) , in order to obtain the BCS formulas for 
the ground state energy and occupation number of levels. We shall follow closely this reference. This continuum limit 
is defined by taking the number of levels N going to infinity while keeping fixed the following quantities: 

g = GN, M/N. (22) 

At the same time the pair energy levels £i will be equivalent to a negative charge density — p(e) located on an interval 
51 of the real axis. The total charge of this interval is given by 

p{e)de = ~. (23) 

We shall suppose that Q is contained inside an interval (— u>,u>). For example, in the BCS model with equally spaced 
levels, oj/2 is equal to the Debye energy. 

The basic asumption made by Gaudin, which is supported by numerical results, is that the solutions of eqs. @ 
organize themselves into arcs (k = 1, . . . , K) , which are piece- wise differentiable and symmetric under reflection 
on the real axis. We shall call V the union of all these arcs, and r(£) the linear charge density of roots in the 
complex plane. Hence, the total number of pairs, M, and the total energy, E, are given by 



r(£) | del = M, (24) 

r 

Jjr(Om = E. (25) 



The continuum limit of eqs. m) is 



*>* P i mm * „, ser , (26) 



e-( Jr ('-( 2G 

which implies that the total electric field on every point of the arcs Tk is null. On the other hand, the values of the 
conserved quantities are given in the continuum by 

The formal solution of eqs. ( p6| ) can be found as follows. First of all, let us orient each arc from the point to 
the point bk, and call Lk an anticlockwise path encircling Tk- We look for an analytic field h(£) outside T and the set 
O, such that 

r(0|dCl=2^(M0-^-(0)de. Cer, (28) 

where h+(£,) and h-(£) denote the limit values of h(£) to the right and left of V. This can be understood using the 
electrostatic equivalence, considering that the electric field presents a discontinuity proportional to the superficial 
density of charge when crossing such surface. Next, we define a function with cuts along the curves r^, by the 

equation 



K 



,1/2 



n(£-°*x£-&*) 



,fe=i 



(29) 



and look for a solution which vanishes at the boundary points of T, in the form 
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HO = R(0 



tp(e) de 



in e-C 

This field has to be constant at infinity, hence the first K — 2 moments of tp must vanish, i.e. 



(30) 



e l <p(e) = 0, 



<£ <K- 1. 



(31) 



This formula is empty if K = 1. The contour integral surrounding the charge density in (pq) can be expressed as 



KO \d?\ 



2?ri 



/z + (Q-fe-(Q 



£ecr, 



(32) 



where Cr is the region outside the curves in T. Using eqs. (129) and (BG), one finds for the principal value of (B2| 



■I 



r{?) \d?\ 1 f dCR(e) f <f(e) de 



2ni 



y(e)fl(e: 



cfe 



/ £ *-V(e)<fe, £ e r, 

JO 



(33) 



where we have deformed the contour of integration L into two contours, one encircling the interval £1 (first term) and 
another one around the infinity (second term). We are assuming that T, and consequently L, do not cut the interval 
fl. The case where L intersects Q actually arises in the equally spaced model, and will be discussed later on. 
Plugging eq. (|33]) into (26), we see that a solution is obtained provided the following equations are satisfied: 



_ g(g) 
R(sY 

K-l P( £ ) 



de 



1 



which have to be suplemented by eq. (|3~i"|), which can be rewritten as 



(34) 
(35) 



E *?&de = 0, 0<£<K-1. 



(36) 



The field h(£) gives also the density of charges r(£), 



and its value is given by replacing (134) into (BOJ), i.e 



r(0 = -K0h Cer, 

7T 



p(e) cfe 



Jo 



(37) 



(38) 



whose value at infinity is —1/2G. Moreover, using eq. (|32|) and performing the same contour deformations that lead 
to eq.(33), one can show that the conserved quantities ( p7|) are given by 



He) 



GR(e) 



p(e') de' 



R(e') e'-e 



G 



(h(s + i0)+h(e-i0)), 



(39) 



so that A(e) is the principal part of h(e) on the set f2. Finally the equations fixing the arcs r& arc the cquipotcntial 
curves of the total distribution, i.e. 



n [ h(0 d? = o, £ e r fc . 

J a k 



(40) 
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VI) The BCS equations 

In this section we shall show how the BCS equations describing the ground state of the model can be derived from 
the formalism of section V jlO| . 

The basic assumption is that for the ground state all the roots En form a single arc, i.e. K = 1, leaving only two 
complex parameters a = a\ and b = b\, which shall be denoted as 

a = s -iA, b = e a + iA, (41) 

where Eq and A > are real parameters. 

From cq. ( |38|) the electrostatic field is given by 

M£) = [« - ^o) 2 + A 2 ]* / p & JL. . (42) 

J si [(s-e ) 2 + A 2 ]3 e-£ 

And the conserved quantities A(e) can be derived from the eq. (|39|). Eq. (|35|) yields the BCS gap equation 

P{£)d£ = -1 . (43) 
Similarly, eq. (E4) becomes the chemical potential equation 



27ri Jl Jq \ 



e - £ 



v /(e-eo) 2 + A2 

while eq. (p5|) gives the BCS expression for the ground state energy, 



p(e)de, (44) 



E = ^- f £ h (t)dt=~ + / (l- £ g |p( £ ) £ &. (45) 

2tt« Jl 4G 7 ^ ^(e - e ) 2 + A 2 J 

Gaudin's paper |l(| contains a misprint in this equation since the term — is quoted as — |g. 

Comparing these equations with the corresponding ones in the BCS theory, we deduce the following relations 
between A, Eq, and p(e), and Abcs (BCS gap), p, (chemical potential), and n(e) (single particle energy density): 



A = 2A BCS , E = 2p, p(e) = -n(^y (46) 

The factor 1/4, in the last term of (|4^), emerges from the normalization 1/2 of the charge density p(e), times another 
factor 1/2, due to the fact that the separation between energy pairs is twice the separation between single particle 
energy levels. 

To determine the equation of the curve T we define the complex variable 



z = V(£-£o) 2 + A 2 = x + iy . (47) 

Using eqs. (Eoh and (E2|) we deduce for T, 



/ de M J, ,2 + A 2 + f ^ 72 • (48) 

Jsi I v / (e-£o) 2 + A 2 4 / + ^ (£ _ £o)2 + A 2A , ,.2 I 



2 

2^ My? 

2 



Eqs. (H^), (fHj), and (EB|) have been derived under the condition that T does not cut the set fl, which requires that 
( [48| ) has no solutions with £ (E Q. In the equally spaced model the curve T may intersect fl, so this case has to be 
treated with some care (see below). 
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Finally, the mean occupation of the energy level Ej can be computed from the formula (Nj) = dE/dsj. In the large 
TV limit (Nj) can be shown to be given by [||, |l(| 

^'K'- tfot )' e£Si - ,4S) 

which agrees with the BCS result Next we shall apply the previous formalism to two models. 

The two-level model 

Let us suppose that we have two energy levels ±£i, which can accomodate at most N/2 pairs each. The charge 
density is given by 

P (e) = j[S(e + e 1 ) + S(e-e 1 )}. (50) 
From eq. (E3) we deduce that the chemical potential vanishes, i.e. Eq = 0, while the gap eq. (Ha) gives 

e\ + A 2 = g\ g = GN. (51) 
This equation implies that g must be greater than E\. The total energy of this solution follows from ([45|), 

E=-^(e 2 1+ g 2 ). (52) 
Integrating eq. p3) we obtain the electric field, and accordingly the value of the conserved quantities 



M0 = ^^P# — A(- £l ) = ^(2 £ ?+ 5 2 ). (53) 
2.9 sf - ^ 8ge! 

Similarly one finds \(si) = — A(— £i), in agreement with the condition = derived in sect. III. The equation of 

the curve (fl8]) becomes in this case 

-+log|^±4 = 0, (54) 

g {x + g) +y 

or equivalently 

* 2+ ^ 2 =thSfer ^ 

where x and y are given by eq. ([!?]) . In Gaudin's paper there is a misprint in eqs. (^) and (|5^). 

From the numerical results shown in fig. |l|, when g > e\ T is an open arc between the imaginary points ±iA in the 
complex £-plane. However when g — E\, in the limit when N — > oo, the arc closes at the origin, and for g < E\ it 
remains closed surrounding the point —£\. 

To treat this case Gaudin considers a general situation where T is a closed curve surrounding a piece uj of f2 and 
defines a function s(£) by the equation |1(J 



Then eq. (^6|) is solved by 



(56) 
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*(0 



1 f p(e) de f p(e) de 
2G 



(57) 



which can be verified using P J T d£' / (£' — £) = in, when £ £ T. For the two-level model s(£) is given by 



and eqs. (53) and <M) yield 



5(0 



1 

2G 



N e x 



2e\ 



(58) 



N 
~2 



M — [ d£ s(0 = / de 2 j o(e) 



(59) 
(60) 



Hence the state for g < ei is a Fermi sea to leading order in N. 

The conserved quantities are derived from an equation similar to (|39|), where h(£) is replaced by s(£) ; i.e. 



7rG 

A(e) = -i— (s(e + iO) + s(e - iO)) , e e w 



(61) 



verifing, as in the previous case, the symmetry A(ei) = — A(— £i). 

To find the curve T we integrate (|5q), and impose the result to be real accordingly to ((56|). Impossing the imaginary 
part of eq. (B6[) to vanish we obtain a diferential equation for the curve T, which is a check of consistency. 



(62) 



n ( — + \ ^ 7TT ) = °' ^ r < 

where Xq is an integration constant, which has a non vanishing value, contrary to Gaudin's assumption. Eq. ( |62] ) can 
be written as 



x 2 + y 2 + e\ = 



2xe\ 



th{2{x-x )/g)' 



£ = x + iy. 



(63) 



To find Xq we have first to look for the point So where T intersects the real axis. From fig. [I], the density of roots at 
this point vanishes, i.e. s(eq) — 0, which implies 



£o = -£i >/l - g/ei . 
The value of x$ is found by imposing that £q belongs to the curve (|63|), namely 



(64) 



9 , £i+£a 

xo = e Q ~ - log 

2 ei - £q 



(65) 



As g approaches £\ from below, both £o and xq go to zero, and the curve ( p3[ ) coincides with the curve ( |55| ) at the 
"critical point" g = e\. 
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Equally spaced model 

Let us consider now a model with N pair energy levels that are uniformly distributed in the interval [— with 
a density p(s) = po 

The distance between the single particle energy levels, d, is given by d = ui/N — 2u>d/N, where cj.d is the Debye 
energy. This is the model used in most of the studies on ultrasmall superconducting grains |35[ ). At half-filling the 
number of pairs is M = N/2, which implies that Eq = (see eq. (ff4|)). On the other hand, the gap eq. ( |43| ) yields the 
standard result H 

w GN 

A= . , 3= , (67) 

smn(l/j) u> 

where g is the usual dimensionless BCS coupling constant. Recall that A is twice the BCS gap Abcs- The total 
energy of the BCS ground state can be derived from (f45h, 



A^ 2 



E =- 4^ V 1+ UJ ' (68) 

while the energy of the Fermi sea state is given by 

^ 2 GN 

EFS = --d~ — - (69) 
The condensation energy is defined as Ec — E — Eps, which in the weak coupling limit, where A <C ui, behaves as 

(™) 

the well known result. 

Next, we shall study the shape of the arcs in this model. The electrostatic field h(£) is given by 



m = « io g «-")^ , -{"+^g , HA» + M ») 

(f + o.)(A 2 + ?w + \J (A 2 + f)(A 2 + u.' 2 ) 



To simplify matters, Gaudin considered the limit A«u, where eq. (|7l|) becomes 



The conserved quantities A(e) are given by 



h(0 = po log L-Vg!±^, |£| « u. (72) 



Mc) = g lo . Kg - ^)(A 2 - eu + y (A 2 + £ 2 )(A 2 + ^ 2 )| J , . 

4 8 |( E + W )(A2 + eu; + ^(A 2 + e 2 )(A 2 + w 2 )| ~ ' 4 g \e + VA 2 +e 2 \ ' 

In the low energy approximation the equation of the equipotential curves passing through the points £o ± iA is 
jiven by 



n 



z + v/a 2 + e 



(74) 
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This curve cuts the real axis inside the interval [— u>,ui], which contradicts the assumptions made so far. From 
numerical studies we can see that for weak couplings only a fraction of energies En form complex conjugated pairs, 
which in the limit N — > oo organize themselves into arcs, while the other energies remain real and located between 
the lower energy levels (see the appendix for more details on this issue). Taking this into account and using contour 
arguments, it can be shown that the BCS eqs. Q), and (Ea) also hold in cases where T cuts the energy interval 

n. 

In the general case the curve T is given by the equation 



K 



. l + £ 2 /A 2 . A 2 
it; arcsm. + w ArgsmhW — g- 



1 -e/u 2 



0. 



(75) 



In the limit A<w eq. (f75j) turns into ([74|) . We show in the appendix that for A > to the arc T does not touch the 
interval [— u>,uj\ and all the roots are complex. This happens for g > go = 1.13459 (sinhl/go = 1). 

VII) Numerical results 

In this section we present the comparison between the analytic results derived in section VI, and the numerical 
ones obtained by solving directly the Richardson's eqs. (j|) for the models considered above. 

Two-level model 

In this model there are two energy levels, ±£1, with a degeneracy N/2. Hence the Richardson eqs. (||) become 



1 

G 



N 

~2 v .-, 



E v 



E v 



M 

E 



En, — En 



(76) 



with M = N/2. The total number of solutions of this system of equations is given by the combinatorial number C^. 
The solution which corresponds to the ground state of the BCS Hamiltonian is the one for which E^G) — > — £1, V/i 
in the limit where G — > 0. For any non zero value of G, all the roots E^ form complex conjugate pairs surrounding 
the lower energy —£\. 

The reduced coupling constant for this case is taken as g = GN . Our numerical results are presented in units of 
£1, which is equivalent to set E\ = 1. 

In fig. H we plot the distribution of the M — 100 roots E^, for three different values of the coupling constant g — GN, 
together with the analytic curves derived in section VI. There is an optimal agreement between the numerical and 
analytic results in the three regimes: normal (g < E\), critical [g — £1), and superconducting (g > £1). As we 
increase the value of M up to 800 pairs, the fit of the numerical data to the analytical curves improves, and in the 
limit M — > 00 the discrete roots collapse into the continuous curves. In particular, the complex roots E max (M) and 
E^ nax (M), whose real part lies closer to £q, approach their continuum limit value, namely 



lim E„ 

M— too 



S (M) =£ + iA, 



where 



so 



A 



g > £l 

-e x y/l - g/ei g<£i 



g > £\ 
g < £1 



(77) 

(78) 
(79) 



in agreement with eqs. (^l|) and (|64|). 

Fig. H shows the following scaling behaviour E max (M): 



E max {M) - (£ + »A) 



(80) 
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FIG. 1: Plot in the £-plane of the solutions of eq. ( |76[ ) for M=100 and three values of g. The continuous lines are the theoretical 
curves (p5J) (with the change of variables fcj)) for g = 1.5, 1; and for g = 0.5. All the energies are in units of e%. 
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FIG. 2: Plots of A — Im E max (M) and eo — Re E max (M) versus M 8 , for the two-level model, eo and A are given in eqs. (178 
and (|79[), while numerical values of 9 are given in table 1. 
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where the exponent 9 depends on the regime of the system. The numerical and analytic results of are given in 
table 1: 



e 


9 = 1-5 £i 


g = 1.0 ei 


g = 0.5 ei 


analytical 
numerical (A) 
numerical (eo) 


2/3 
0.635 
0.672 


1/3 
0.329 
0.328 


1/2 
0.446 
0.513 



Table 1.- Numerical and analytical values of the exponent 9 appearing in eq. (p0[). The numerical values are 

extracted from the real and imaginary parts of E max (M). 

The analytical value of 9 can be derived as follows. For M sufficiently large we can write 

E max (M) = {e +iA) + 5Z, (81) 
where <5£ is a small number. Using eqs. (p3|), ([56]) and (p8|), the density of roots r(E max ) behaves as 

f 1/2 g>ei 

r(E max )~AM(60", v={ 2 g = £l , (82) 

[ 1 .9 < £i 

where A is a parameter which only depends on g and ei. The magnitude of St; can be fixed by imposing the existence 
of at least one root in the interval (eq + iA, Eq + iA + 6£), namely 



1 = r(E max ) \5£\ ~ AM \5tf +1 — > <5£ L- , (83) 

which is the desired result fl80| ) with 9 = 1/(1 + v). 

As we have seen in fig. [l], the roots E^ lie on the analytic curves to a very good degree of approximation. A further 
check of the analytic solution is to compare the density of roots r(£), which is given by r(E) = \h(E)\/-K for g > E\ 



(see eq. (p7|)), and by r(£) = |s(£)| f° r 3 < £ i ( see e Q- (56)). In fig. |3j we plot the fraction of roots per unit length 



r(E)/M as a function of the arc length s£ [—1,1], where the s = ±1 corresponds to the end points of the arc T. The 
agreement is extremely accurate for M = 800. 

The value of the conserved quantities obtained in eq. (|5^) also agree with the numerical ones in the large N limit, 
in the three different regimes. 

Finally, it is interesting to study the behaviour of the energy as a funtion of the number of pairs M. The leading 
term is given by eq. Js^) for g > E\ and by eq. (p0|) for g < E\. In fig. || we plot eoo — e(M) versus 1/M, where 
e(M) — E(M)/M is the energy per pair. The linear behaviour of the numerical data means that the next leading 
correction to E(M) is a constant term, E\, whose value is given by the analytic expression 



-^/M = e.-e(M) = {fe-^)/« «|*. (84, 

where A and Eq are given by eqs. ( |5l| ) and (|64] ) respectively. The formula for g > E\ has been derived by Richardson 
in ref. g, studying the next to leading order of the electrostatic model. The formula for g < e% is an educated guess 
which needs a proof. At g = £\, the numerical data show a large deviation, for small M, from the formula (|84|) , 
namely £\/M . There are theoretical reasons to believe that the scaling behaviour of E\jM at g = E\ , is given by the 
following formula (see fig. ^): 



-E 1 /M = e 1 (±-j^,Y A = 0.62258, (85) 

which would imply that the next-to-next leading correction of the total energy goes as M -1 / 3 , rather than 
M~ x . In fact, the formula given by Richardson B for the 0(M~ 1 ) energy correction, which is given by 
E-2 = —g(3g — A)(g — A) / (4N A 2 ) , breaks down at g = E\ where the gap vanishes. The derivation of eq. ( p5[ ) shall 
be given elsewhere. 
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FIG. 3: Density of roots r(E)/M as a function of the arc length. The continuous curves give the analytic result while the 
step like curves are the numerical results obtained for M= 800 and three values of g. The total length of the curve has been 
divided into 50 parts. The height of each part is the fraction of roots £" M per unit length SE. 
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1/M 



FIG. 4: Plot of eoo — e(M) versus 1/M for the two-level-model and different values of g. The continuous lines are given by 
eqs. (Q) for g / ei, and eq. ( pBj ) for g = e\. 
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FIG. 5: Evolution of the real and imaginary parts of E )1 (g), in units of d — uj/N, for the equally spaced model with M — N/2 — 8, 
as a function of the coupling constant g. For convenience the energy levels are choosen in this figure as ej = 2j. This figure 
confirms numerically the behaviour found in reference H. 



Equally spaced model 

This model is defined by N non degenerate energy levels Sj = d(2j — N — 1), j — 1, . . . , N, where d = u/N is the 
single particle level spacing. The first numerical study of this model, using the exact solution, was done by Richardson 
in reference Q, where he introduced new set of variables, related to the roots E^, in order to handle the singularities 
arising at particular values of the coupling constant g (see below). This work was revised recently in connection to 
ultrasmall metallic grains in |56| 

The ground state of the system, at half-filling M = N/2, is the one for which lim g ^o E^(g) — ► (fi = 1, . . . , N/2). 

The reduced coupling constant is given by g = G/d = GN/uo in this section. And all our numerical results are 
presented in units of u, which is equivalent to set lo = 1. 

In fig. H we plot the real and imaginary parts of E^, for a system with M = N/2 = 8 pairs. For small values of 
<7, all the pair energies are real and below the corresponding energy levels at g = 0. At a certain value g\, the two 
closest pairs to the Fermi energy, namely Em and Em-i, coincide at Em-x and then, for larger values of g > g\ 
they form a complex conjugate pair, while the rest of the pair energies remain real. At higher critical values of g the 
same phenomena happens for the other energies, until all of them become complex. At the critical values of g the 
Richardson eqs. are singular, but the singularities can be resolved by the change of variables introduced in ||. 

In fig. ^ we show the distribution of roots for M = 100 pairs, together with the theoretical curves obtained in 
section VI, for three different values of g. For g = 0.5 and 1 there is a fraction of roots which become complex, falling 
into the arcs T described by eq. (|7^) , while the real roots extend from the cutting point of the arc T with the real axis 
down to the energy — u). For g = 1.5 all the roots are complex and fall into the theoretical curve (|7^). The numerical 
data are consistent with the value go — 1.13459 above which all the energies become complex. 

As in the two-level case, the root E max (M) closest to eo + iA satisfy the scaling law ( |8C| ) with an exponent 9 = 2/3. 
This result can be proved along the same lines as for the two-level case. In table 2 we give the comparison between 
the theoretical and numerical values of 0, and in fig. |l^ we show the scaling behaviour for three values of g. 
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FIG. 6: Plot of the roots for the equally spaced model in the complex £-plane. The discrete symbols denote the numerical 
values for M = 100. The continuous lines are given by eq. ([T5|). All the energies are in units of w. 
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Arc Length 



FIG. 8: Density profile of the complex roots for the equally spaced model with M — 800. Notations are the same as in fig. 



6 


.g=1.5 


.9=1 


g = 0.5 


analytical 
numerical (A) 
numerical (eq) 


2/3 
0.646 
0.658 


2/3 
0.645 
0.659 


2/3 
0.642 
0.661 



Table 2.- Numerical and analytical values of the exponent 8 appearing in eq. 

notation is the same as in table 1. 



(gCl) for the equally spaced model. The 



In fig. U we give the fraction of roots per unit length r(E) /M for M = 800 and three values of g. Again, the 
numerical data fit accurately the theoretical result r(E)/M = \h(E)\/w, where h(E) is given by eq. (|7l|). 

A further check is presented in fig. ^, where the distribution of the conserved quantities \{£i) for M — 20 is 
presented, together with their continuous limit curves given in eq. (|73|). A very good agreement is achieved even for 
such a small number of pairs. 

We have finally studied the condensation energy Eg/M as a function of M, obtaining the behaviour 
Eq{M) = eooM + Ei, where E\(g) is plotted in fig. [h]. An analytic expression for Ex(cj) is not known In 
order to understand the numerical results we have made two fits. The first one is based on the result obtained for the 
two-level case, namely 



- E x {g) = gu - - , (86) 

where A is given by eq. (|67]). Fig. |l^ shows that (|8^) is a good fit for large values of g, where all the roots are complex, 
forming a single arc. For g < go — 1.13459 there is a fraction M comp i ex /M of complex roots, and a fraction M rea i/M 
of real roots. This suggests the following mixed fit: 



-Mg),.= (9-^)^^ + (A + Bg)g^f, 

(87) 

A = 0.6735, B = -0.1729, 

which describes pretty well Ei(g) in the whole range of g's. In the weak coupling regime there have been proposed 
several formulas for the finite size corrections of the condensation energy Eq using the DMRG method (33]] and the 




FIG. 10: Plot of —Ei(g) versus g for the equally spaced level model. The fits of the numerical data are given by eqs. (g6|) and 
(B7p in units of ui. 
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exact solution j34j. In the latter reference it is shown that for couplings g > g* = 1/lnN, the condensation energy 
Ec can be splitted in two contributions, one due to the complex energy pairs ( the BCS term) and another one due 
to the real roots ( perturbative term). It would be interesting to clarify the relation between formula ( |87| ) and eq.(24) 
in reference 134] . 



VIII) Conclusions and prospects 

In this paper we have reviewed and completed Gaudin's continuum limit of the exactly solvable BCS model, and 
compared its predictions with the numerical solutions describing the ground state of the two-level and the equally 
spaced models. We have confirmed previously known results |(| |25|, |3(| p7| , and obtained new ones, which we list 
below: 

• Analytic and numerical determination of the curves T, where the roots of the Richardson equations lie, 
together with their density along T. For the equally spaced model the curves T are given for arbitrary values 
of the coupling constant g, and for the two-level model in the normal regime, i.e. g < ex, the analytic curves T 
differ from those in jio) . 

• Analytic expression of the eigenvalues Xj of the conserved quantities found by Cambiaggio, Rivas, and Saraceno 
[ [TTf , which fit accurately the numerical data, even for a few levels. 

• Study of the scaling properties of the ground state energy E(M) with the number of pairs M which, in the large 
M limit, behaves generically as E(M) = e^M + Ei +0(1/M). For the two-level model, in the superconducting 
regime, we find an agreement with the analytic result for E\ [||, while in the normal regime we have guessed an 
analytic expression for E\ , which fits perfectly the numerical data. In the critical case there are large deviations 
in E{M), which can be explained by a term (^(M" 1 / 3 ). In the equally spaced model the condensation energy 
also receives a constant contribution, Ex, which is fitted with a "phenomenological" formula which combines 
the effects of the complex and real roots. 

• Scaling behaviour of the roots E max (M) which lie closest to the end points of the curves T. This scaling can be 
characterized by a critical exponent 9, whose analytic values fit well the numerical ones. In the superconducting 
regimes (i.e. Vg > for the equally spaced-model and g > e± for the two-level model) we find 9 = 2/3, while 
in the normal (critical) regimes of the two-level model we find 9 = 1/2 (1/3). In this manner the exponent 9 
characterizes the nature of the ground state. An interesting question is wether 9 shows up in physical observables. 

• Generalization of Gaudin's conformal mapping for the equally spaced model using elliptic functions, which 
degenerate into trigonometric ones in the limit A <C to. This suggest the existence of a rich mathematical 
structure. 

Finally, we shall mention some problems which are worth to investigate along the lines suggested by the present work: 

• Study of the excited states and the finite size effects in a more geometrical way, completing the work initiated 
by Richardson in ^ . 

• Quantum dots at finite temperature is another difficult problem, which has been treated combining several 
techniques depending on the temperature regime |26[ p7| , |38| . Is it possible to find a thermodynamic Bethe 
ansatz for this system? 

• Study of solutions of Richardson equations with several arcs, i.e. K > 1 in the notation of section V. For the 
standard BCS model they must describe very high excited states formed by separate condensates in interaction. 
This case may be relevant to systems such as arrays of superconducting grains or quantum dots. From a 
mathematical point of view, the cases with K > 1 seem to be related to the theory of hyperelliptic curves and 
higher genus Riemann surfaces, which may shed some light on this physical problem. 

• Generalization of Gaudin's equations to Richardson models based on arbitrary Lie groups Q [|l6| . The standard 
case corresponds to the choice Q — SU(2). The electrostatic model for generic Q has charges labelled by vectors 
of dimension rank Q. It is naturally to expect that the roots E£ (a — 1, ... , rank Q) should fall into several arcs 
labelled by the index a — 1, . . . , rank Q. 
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APPENDIX: CONFORMAL MAPPING FOR THE EQUALLY SPACED SYSTEM 



In this appendix we study the analytic structure underlying the shapes of the curves defined by eqs. (74) and (ul 



i) Case A < u> 

In this limit it is convenient to make a conformal mapping jl^ from the upper half-plane Im£ > 0, with a cut in 
the segment (0, iA), into the half-strip A > 0, — n/2 < p < tt/2, 

£ = iA cosh(A + iji). 



In the plane u = A + ip, the curve (74) becomes 

A tanhA + p cot p = 0. (89) 

Figures [ll| and [l^ show, in the u and £ planes, the arcs defined by eqs. ( p9| ) and (frij). The point O (£ = e + iA) is 
mapped into the origin (it = Q), where the field /i vanishes. In fact h(£), as a function of it, is simply h = —2pou. The 
point O' (£ = £o — *A) is the mirror reflection of O. The arc OA is the one that corresponds to the actual numerical 
solution for Im£ > 0, while its mirror image, O'A, is the one for Im£ < 0. We thus see that the whole arc OAO' cuts 
the real axis at the point A with an energy 

£a = £a — A sinhAo, Ao tanhAo = 1. (90) 

These results are expected from numerical studies, which show that for weak couplings only a fraction of energies E^ 
form complex conjugated pairs, which in the limit N — > oo organize themselves into arcs, while the other energies 
remain real and located between the lower energy levels ( see fig. ||). We can check this result computing the number 
of real and complex roots E^. From eq. ( p0[ ) we deduce that the real roots occupy the interval [— u, — AsinhAo], with 
a density which is 2po- Hence the total number, M rca i of real roots is 

M real = 2p (u> - A sinhA ). (91) 

Since the total number of roots is M = M rea j + A/ comp i ox = N/2 = 2pou>, we deduce that 

M comp iex = 2Ap sinhA . (92) 

This result agrees with the one obtained integrating the field h(£) along the curve L surrounding T (recall eq. fl44|)). 
This integration is more easily done in the it-strip 

M comp i ox = / fe(f)d£ = 4x / {-2p u)(iAsmhudu), (93) 

JL JO 

where the factor 4 comes from the contribution of the up and down, and left and right pieces of the contour L. Using 
contour arguments, it can be shown that the BCS eqs. ((43|), (|44|), and ( {45| ) also hold in cases where T cuts the energy 
interval Q. As a side comment, we observe that Mcompiex = ^ sinhAo, agrees with the heuristic argument stating that 
there are roughly Ascs/d energy levels around the Fermi level that are strongly affected by the pairing interactions 
f39f . See also reference |54| for an approximate equation giving the number of complex pairs ( eq. (49)). 



ii) Generic case 

In order to compare with the numerical results presented in the section VII, it is convenient not to make the 
approximation A -C ui. The appropiate change of variables that generalizes (pq) is given by 

iA cosh it 

£ = == , (94) 



— sinh it 

to 
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FIG. 11: Half-strip in the u — X + ifi plane that is conformally equivalent to the upper half £-plane, with a cut in the interval 
[0, Eo + iA]. The arcs OA and OB are given by eq. (^). This figure is borrowed from reference JlQ ], 
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FIG. 12: Location of the roots in the £-plane. The complex roots condense in the arc OAO', given by eq. (^), while 
the real ones are distributed in the segment X' A, with a real root between every two energy levels. The point A is given by 
eq. (pel). This figure is borrowed from reference 
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in terms of which h — —2p^u. The integral J h(£)d£ can be more easily done in the u plane, yielding for the equation 
of the curve T 



i u cosh u 



A 



sinh u 



!W . / A 

= — arcsin — sinh u 

2 A \lo 



= 0, 



(95) 



which in the £-plane gives eq. (|75|), namely 

K 



/l + £ 2 /A 2 /A 2 + C 2 

it; arcsin W — — - + to Argsinh- 



i -e/u 2 



0. 



(96) 



In the limit A <C to this equation turns into eq. (|7J). We can ask when this curve cuts the interval [— w,u)]. A 
convenient parametrization of the cutting point is (recall eq. (|90|)) 



u = - 



A sinh Aq 



, A 

1 + — cosh Ao 



(97) 



Plugging < p7\) into ( |9 6| ) we deduce the equation for A , 

Aq sinh A 



I ( — coshAo 



uo . /A 
= = — Argsinh — cosh Ao 
2 A \ UJ 



which only has real solutions provided A < oj. When A > uj the arc V does not touch the interval [— uj,lo] and all 
the roots are complex. This happens for g > go — 1.13459. This value of go is obtained solving the equation 
sinhl/go = 1. 



To complete this appendix we shall show that eq. (94) yields a conformal mapping that generalizes the one found 
by Gaudin for generic values of A/a;. To do so, let us first perform the change of variables 



u = i 



and introduce the parameters 



in terms of which ( |94| ) becomes 



k = 



VA 2 + w 2 



k' = 



VA 2 + uj 2 



= Vl-fc 2 



a/I - k 2 sin 2 (j> 

Using the elliptic functions with modulus k, and conjugate modulus fc', we define a new variable z through 

sin0 = sn(z,fe), 
cos0 = cn(z,/c), 



(99) 



(100) 



(101) 



(102) 



'l — k 2 sin 2 cj) = dn(z,fc). 
The relation between z and <p is given by the equation 



z{k,<f>) 



r-q> 

Jo \fl 



\]\ — k 2 sin 2 %j) 



(103) 
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while the inverse function is denoted as <fi = & m (z, k). The relation between £ and z reads 

., sn(z,k) 



£ = iAk 



dn(z, k) 



(104) 



Fin ally, using the relation sn(K + v,k)/ dn(A" + v,k) = cn(v, k)/k', where K is the half-period of the elliptic integrals, 
eq. (104) turns into 



£ = iA cn(v, k), 

where z = K + v. The relation between the variables u and v is given by 

2q n 



IT 

-— + an i 



Sill 



mr(K + u) 
2~K 



(105) 



(106) 



whe re q — e ' K . In the limit A/u> « 1 we get k — > 0, K — ► 7r/2, and g — > 0, which implies that iu —> v and 
eq. ( |i05l) re duces to ©. 

Eq. (105) gives a conformal map of the rectangle —X < Rc v < K, < Im v < K' onto the upper half-plane 
Im £ > 0, with a cut along the linear segment < Im £ < A (see reference jy] and fig. [l3|). Some examples of this 
mapping are given by 



v = 0, ±K, ±K + iK' 



£ = zA,0,±c^, 



(107) 



which imply the following identifications between the sides of the rectangle in the u-plane and the relevant energy 
regions in the £-plane: 



cut 


Re£ = 
< Im £ < A 


-K < Re v < K 
Im v = 


lower band 


-u < Re £ < 
Im£ = 


Re v = -K 
< Im v < K' 


upper band 


< Re £ < lu 
Im£ = 


Rev = K 
< Im v < K' 


off band 


|Re || > u 
Im£ = 


—K <Rev<K 
Jmv = K' 



Table 3.- Correspondences of the boundary regions of the conformal map (105). 



